function traces = add_response(traces)%,... traces:irisFetch trace struct
    %water_level) %max magnification relative to min magnification
%convolve instrument response
% Josh Crozier 2019 (some snippets from various GISMO functions)

for tracein = 1:length(traces)
    %make even
    if mod(traces(tracein).sampleCount,2) == 1
        traces(tracein).data = traces(tracein).data(1:end-1);
        traces(tracein).sampleCount = length(traces(tracein).data);
        traces(tracein).endTime = traces(tracein).startTime...
            + traces(tracein).sampleCount/traces(tracein).sampleRate;
    end
    
    %some seem to have different names for terms
    if isfield(traces(tracein).sacpz, 'constant')
        constant = traces(tracein).sacpz.constant;
    else
        constant = traces(tracein).sacpz.k;
    end
    if isfield(traces(tracein).sacpz, 'zeros')
        zero = traces(tracein).sacpz.zeros;
    else
        zero = traces(tracein).sacpz.z;
    end
    if isfield(traces(tracein).sacpz, 'poles')
        pole = traces(tracein).sacpz.poles;
    else
        pole = traces(tracein).sacpz.p;
    end
        
    
    % CALCULATE COMPLEX RESPONSE AT SPECIFIED FREQUENCIES
    ws = 2*pi*(0:length(traces(tracein).data)/2).'...
        *(1/length(traces(tracein).data)/2)...
        *(traces(tracein).sampleRate/2);
    resp = freqs(constant...
        *poly(zero),...
        poly(pole), ws);
    %convert to vel response
    resp = resp./(1i*ws);  %no 2pi to match iris/obspy
    %make 2-sided response
    resp(1) = 0;
    %ws = [-ws(length(traces(tracein).data)/2 + 1:-1:2);...
    %    ws(1:length(traces(tracein).data)/2)];
    resp = ([real(resp(end:-1:2)); real(resp(1:end-1))]...
        + 1i*[-imag(resp(end:-1:2)); imag(resp(1:end-1))]);
    
    %water level
%     abs_resp = abs(resp);
%     max_resp = max(abs_resp(abs_resp>0));
%     resp(abs_resp < max_resp/water_level) = ...
%         resp(abs_resp < max_resp/water_level)...
%         *max_resp/water_level...
%         ./abs_resp(abs_resp < max_resp/water_level);

    
    %debug plots
    %{
%     figure(1)
%     subplot(2,1,1),loglog(1./ws(ws>0),abs(resp(ws>0)))
%     subplot(2,1,2),semilogx(1./ws(ws>0),angle(resp(ws>0)))
    
    figure(2)
    subplot(2,1,1),plot(traces(tracein).data)
    %}

%     window = tukeywin(traces(tracein).sampleCount,0.05);
%     traces(tracein).data = detrend(traces(tracein).data).*window;
    traces(tracein).data = ...
        real(ifft(ifftshift(fftshift(fft(traces(tracein).data)).*(resp))));

    %debug plots
    %{
    subplot(2,1,2),plot(traces(tracein).data)
    xlabel('x')
    %}
end

end

